% MOC DEMO TO DO:
%  1) find a way to approximately verify the code (w/ MCNP or SN)
%  2) precomputation of various factors
%  3) add a cartesian mesh ability for creating regions (w/o circles for
%     now)  This will be very nice for benchmarking, debugging, and
%     visualizing flux maps
%  4) add reflective boundary condition (cyclic tracking verify???)
%  5) add ability to visualize outgoing angular flux distributions at a
%     global edge
%  6) volume-corrected tracks



% main file
% solves for the scalar flux in a 2 region pin cell
% based on code from a. hebert
clear;% clc;

% define the problem
Pitch           = 1.26;
Radii           = [0.54];
NumSpace        = 8;
NumAzimuth      = 8;
NumPolar        = 1;
% region data.  Inner to Outer
% SigmaT          = [ 0.10 0.01 1.00 ];
% SigmaS          = [ 0.03 0.01 0.50 ];
% Q               = [ 1.00 0.00 0.00 ];
SigmaT = [1.0;1.0]; %;1.0];
SigmaS = [0.0;0.0]; %;1.0];
Q      = [1.0;1.0]; %;1.0];
BC              = 0;    % 0=vacuum, 1=reflective (requires cyclic tracking)
MaxIt           = 50;   % Maximum Inner Iterations
EpsPhi          = 1e-4; % relative scalar flux tolerance

% compute the tracks
%  figure(2)
 track   = tracker(Pitch, Radii, NumSpace, NumAzimuth, NumPolar);
% for i = 1:12
%     plotter(track, Pitch, Radii, i)
%     hold on
% end
% solve or the region fluxes
 flux    = mocsolver( track, BC, Q, SigmaT, SigmaS, MaxIt, EpsPhi );
 disp('flux=')
 flux'
 % 7.4613e-01   5.1367e-01   vs  7.2526e-01   5.3059e-01

